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Abstract 

Computing response functions by following the time evolution of superoperators in Liouville 
space (whose vectors are ordinary Hilbert space operators) offers an attractive alternative to the 
diagrammatic perturbative expansion of many-body equilibrium and nonequilibrium Green func- 
tions. The bookkeeping of time ordering is naturally maintained in real (physical) time, allowing 
the formulation of Wick's theorem for superoperators, giving a factorization of higher order re- 
sponse functions in terms of two fundamental Green's functions. Backward propagations and the 
analytic continuations using artificial times (Keldysh loops and Matsubara contours) are avoided. 
A generating functional for nonlinear response functions unifies quantum field theory and the 
classical mode coupling formalism of nonlinear hydrodynamics and may be used for semiclassical 
expansions. Classical response functions are obtained without the explicit computation of stability 
matrices. 
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I. INTRODUCTION 



An important ingredient in many-body theories is the abihty to factorize averages of 
products of a large number of operators into products of averages of pairs. This Wick the- 
orem is common to the broad arsenal of techniques used for the treatment of quantum and 
classical systems alike. Quantized fields are used e.g. in Green function perturbation theory 
of many identical bosons or fermions;^^'^'^'^'^'^ Time Dependent Hartree-Fock (TDHF) and 
Time Dependent Density Functional (TDDFT) equations of motion of many electron sys- 
tems^ and the Hartree-Fock Bogoliubov equations for superconductors and Bose Einstein 
condensates^. Classical fields are considered in mode coupling theories of nonlinear hydro- 
dynamics of fluids and glassesi^^ii; cumulant (1/A^) expansions for short range interactions 
in fluids, and Gaussian models of spin Hamiltonian o^^i^^'^i^^i^^'^i^^i^'^i''^^ . 

Green function perturbation theory forms the basis for the powerful Feynman diagram- 
matic techniques widely used in the description of many-particle systems^i^iMi^iSi^. This 
formalism is based on expressing quantities of interest as time- ordered expansions. Equilib- 
rium and non equilibrium Green function techniquea^i^^ employ various types of contours 
which, in effect, transform the computation to a time ordered form in some artificial (un- 
physical) time variable along the contouri22iSL2^. 

The primary goal of this article is to demonstrate that the description is greatly simpli- 
fied by employing superoperator algebra and computing response functions using the density 
matrix in Liouville space^^'^^'^^^. One of the rewards of working in the higher dimensional 
Liouville space is that we only need consider time ordered quantities in real (physical) time 
and Wick's theorem therefore assumes a particularly compact form; no special contours or 
analytic continuations are necessary. The Hilbert space description requires a sequence of 
forward/backward propagations as opposed to the all-forward representation of response 
functions in Liouville space2LSSi22i^2i^. The superoperator approach provides a unifying 
framework applicable to quantum and classical systems, with and without second quanti- 
zation. It thus connects field theories with classical mode coupling theories of fluctuating 
hydrodynamics. Semiclassical approximations are developed directly for nonlinear response 
functions (i.e., specific combinations of correlation functions) rather than for individual cor- 
relation functions, which do not have a natural classical limit and their semiclassical approx- 
imations are thus ill defined. Recent interest in multidimensional Raman techniques gen- 
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erated considerable activity in modelling multitime correlation function o^^i^^i^^i^^i^^i^^i^^i^^i^° . 
The mode coupling simulation of correlation functions using Langevin equations poses many 
problem ati c. These difficulties disappear by modelling the entire response where the clas- 
sical limit is uniquely and unambiguously recovered. The present formalism shows how 
nonlinear response functions may be expressed in terms of lower order response of collective 
variables^^i^L^i. 

In Section II we discuss two strategies for simulating response functions. The first, based 
on the wavefunction in Hilbert space, does not maintain a full bookkeeping of time ordering 
whereas the second, based on the density matrix in Liouville space doesM'^^^. A detailed 
comparison is made of the physical insight and the numerical effort required in both pictures. 
These results form the basis for developing the many-body Green function perturbation 
theory in Section III. Using a generalized superoperator generating functional, we obtain a 
time ordered perturbation theory of elementary Liouville space operators, and derive Wick's 
theorem for Boson field superoperators in Section IV. These results are used in Section V to 
derive a semiclassical expansion for response functions which in the classical limit recovers 
mode coupling theory. The extension to Fermion fields is made in Section VI and our results 
are summarized and discussed in Section VII. 

Wick's theorem is based on a perturbative expansion around a quadratic Hamiltonian 
and is thus limited to physical situations when this is a good reference for the actual dy- 
namics. It is given for Boson fields in Section IV using a closed expression for a generating 
functional, and for Fermion fields in Section VI. In Section V we explore it in coordinate 
space without using second quantization. Section II introduces the notation and reviews 
previous results. The superoperator algebra of Section III was used earlier for specific ap- 
plications (time dependent Hartree-Fock, fifth Raman spectroscopy).— i^^i^^^ This section 
recasts these earlier results in a more general and compact notation that sets the stage for 
the subsequent sections. 

II. LIOUVILLE VS. HILBERT SPACE DESCRIPTION OF QUANTUM NONLIN- 
EAR RESPONSE 

Partially Time Ordered, Wavefunction Based Expansion of Response Functions 

We consider a material system with Hamiltonian H, coupled to an external driving field 
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E{t) by the interaction 

Hintir) = -AE{t), (1) 

where A is a general dynamical variable. For clarity we assume a scalar field; Extension 
to vector fields is straightforward by introducing tensor notation. The total Hamiltonian 
Ht{t) is given by 

HT{T) = H + H,nt{r). (2) 

We shall be interested in the expectation value of an operator B of the driven system 
at time t. For a system described by a wavefunction | ipj{t)) this is given by S{t) = 
{ip j{t)\B\%l) j{t)) . A perturbative calculation of | then gives to ra'th order in the field 

n 

Sf\t)Y.{^t\i)\B\i^f~"'\t))- (3) 

m=0 

Here iV'j'"'') denotes the wavefunction to m'th order in Hint. If the system is initially in a 
mixed state (e.g. Canonical distribution) where state |j) is occupied with probability Pj, 
we need to average Eq. over that ensemble 

S^-\t) = Y.P'jSf\t). (4) 

Time dependent perturbation theory gives for the linear response^^ 

Sf\t) f dn{^,\U^{t - n)BU{t - n)A\^,)E{T) + c.c. (5) 

Here \ipj) = \'ipj{0)) and U{t) is the retarded evolution operator in Hilbert Space which 
propagates the wavefunction forward in time 

f/(r) = e(r)exp(^-^/7r), (6) 

whereas the advanced Green function 

f/t(r)=e(r)expQ/7r), (7) 

is responsible for backward propagation. 9{t) denotes the Heavyside function (0 for r < 0, 
1 for r > 0). 

For the third order response which describes many of the most common nonlinear spec- 
troscopiesS-, we obtain 

-^f W = (i-)' f r r dr^Ra{t,r^,T2,T^)E{n)E{r2)E{T:,) (8) 

+ (^1 f dn r dT2 r dT3Rb{t,T3,T2,Ti)E{Ti)E{T2)E{T3)+C.C 
J —oo J —oo J—oo 
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where 



i?a(r4,r3,r2,ri) = {iPj\U\r3^)AU\ns)BUin2)AUir2i)A\^j) (9) 

Rcir,,r3,r2,n) = {^j\U\u,)BU{n3)AU{T32)AUir2i)A\^j) . 

and we have defined T4 = t and Tij = ti — tj. These equations represent a time loop of 
forward and backward propagations^^. Eq. may be alternatively recast using correlation 
functions 

Ra{rA,r3,r2,n) = {^,\A{r,)B{u)A{T2)A{n)\^,) (10) 

Rb{n,r3,T2,n) = {^j\A{n)A{T2)B{T4)A{Ti)\^j) 

Rc{U,r3,T2,n) = {^j\B{u)A{T3)A{T2)A{ri)\^j) 



where we denote operators in the Heisenberg picture by ( " ) 

i(r) = U\r)AU{T) (11) 

The time variables of Rc in Eq. (jH)) are fully time ordered (ri < T2 < t-^ < t). However, this 
is not the case for Ra and Rb. By breaking the integrations into various segments we can 
maintain full time ordering, and recast Eq. (jH)) using a response function. This will be done 
next through the density matrix expansion. 

Time-Ordered Expansion: Response Functions 

Rather than using a wavefunction, the state of the system can be described by its density 
matrix, defined as 

Eqs. (121) and ^ can be alternatively recast in the form 

=Tr [5p(")(t)] , (13) 

where 

= E E pM"'\t)){i^t"'\t)\ (14) 

j m=0 
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is the density matrix expanded to the n'th order in Hint- The expectation value of B to n'th 
order in the field is obtained by computing the density matrix to n'th order. This givea^ 



dr, 



n—l 



dri 



R^''\t, r„, r„_i, . . . , n)E{n)E{r2) . . . E{Tn). 



(15) 



Here i?*^"^ is the n'th order response function 

/ „• \ n 

i?^"Hr„+i...ri)= y Tr{[...[[5(r„+0,^(^n)],i(r„-i)]...,i(rOK4- (16) 
which can be alternatively recast as 

/?(")(r„+i . . . n) = Q"Tr{fi(r„+i) [i(r„), . . . , [A{r,), [i(n), Pe,]] •■■]}• (17) 

Note that the time variables Tj in Eq. (jHJ are not time-ordered. In contrast, the complete 
time ordering in Eq. p5|) makes the density matrix description most intuitive and directly 
connected to experimenlj^. 

In the density matrix formulation we maintain a simultaneous bookkeeping of the inter- 
actions with the ket and with the bra. This is why Eq. (fTTj) has 2^ terms, each constituting 
a distinct Liouville space pathway. The wavefunction calculation, in contrast, focuses on 
amplitudes and the various time orderings of the ket and the bra interactions are lumped 
together. Eq. Q has thus only n + 1 terms. The different terms in this case simply reflect 
the order of the interactions within the bra and within the ket (but not the relative time or- 
dering of bra and ket interactions!). When the system interacts with a thermal bath, the 2"' 
terms in Eq. ()17|1 represent very different physical processes and their separate treatment is 
absolutely crucial. The density matrix separates these terms directly and naturally without 
the need for any change of time variables. 

The quantum nonlinear response function R^'^'^ is given by a combination of (n + 1) order 
correlation functions. Response functions provide a natural link between theory and exper- 
iment^, is a purely material quantity which contains all the necessary information for 
describing n'th order response. It is independent of the details of a particular measurement, 
(e.g. temporal sequences of pulses as well as their frequencies and wavevectors) . The fleld 
envelopes enter through the multitime convolutions in Eq. (P3|). When S**-"^ is calculated 



in terms of the wavefunction without using response functions (Eq. we need to repeat 
the calculation for every new realization of the field, i?'-"^ is therefore a compact and eco- 
nomical way for clarifying the fundamental relationships among various techniques and their 
information content. Since the nonlinear response functions are successively probing higher 
order correlation functions, they necessarily carry additional information as the order n is 
increased. 

Forward/Backward vs. All- Forward Representation of Response Functions 
The expression for the response obtained by expanding the density matrix in powers of 
the external field (Eq. fTT)|l . separates naturally into several contributions, each representing 
a distinct time-ordering of the various interactions. The time variables appearing in Eq. (|15p 
are chronologically ordered and represent successive interactions with the field. In contrast, 
the time variables in the wavefunction description are not fully ordered and consequently 
have a much less transparent physical interpretation, i?^") has 2" terms (Liouville space 
pathways) in the density matrix description (Eq. (fTBjl ) but only 2n terms using wavefunctions 
(Eq. (jH))). In practice we need only compute half of the terms (2^"~^^ vs. n) since the other 
terms are their complex conjugates. 
For the linear response Eq. (fTBj) gives 

R^'\r2,n) = l:Y.Pj{^j\U\T2i)BU{T2i)A\ij,) + c.c. (18) 

j 

The third order response is similarly given by 

,T3,r2,ri) = ( -) ^i?f)(r4,r3,r2,ri) + C.C.. (19) 

R?{r4,Ts,T2,n) =J2Pj{^,\U\T2i)AU\rs2)AU\m)BU{ni)A\^j) (20) 

j 

Rf{n,T^,T2,n) = Y,P^{ijj\U{r2i)AU\r,,)AU\n^)BU{r,2)A\i;j) 
j 

R^s\r,,T,,T2,n) = J2PA^j\Uir,MUHr2i)AU^in2)BUiT,,)A\^^) 
j 

R?\r,,T,,T2,n) =J2PA^j\UHr4i)BUir,,)AU{T,2)AUiT2i)A\^^) 
j 

Unlike Eq. (jH)), Eq. (fT3j) allows us to define a response function since it is fully time ordered. 
Note that R4 = Rc, and Ra + Rb correspond to Ri + R2 + R3. 
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Eqs. (ITHj) and (fT^ can be calculated by either expanding the correlation functions in 
eigenstates or using wavepackets in the coordinate representation. Semiclassically it is pos- 
sible to expand | ipjit)) in coherent states | ^Jjit)) = J J dpdq \ pq)(pq | ''Pjif)). Each Rj 
may thus be computed as an average given by a sum over trajectories moving forward and 
backward in time as given by the various U and f/^ factors, respectively. Coherent states 
provide an over complete basis set^. Powerful semiclassical approximations were developed 
for carrying out this propagatio n^^i^^i^^i^°i^^i^^ . 

In Eqs. p9|) and ()20j) we used the density matrix to derive formal expressions for the 
response functions, but for the actual calculation we went back to the wavefunction in 
Hilbert space. Since quantum mechanics is usually described in terms of wavef unctions, 
wavepacket and semiclassical descriptions are normally developed for wavef unctions. It is 
possible however to construct an alternative forward propagating wavepacket picture by 
staying with the density matrix in Liouville space all the way. To that end we represent the 
time dependent density matrix as 



The first equality is the common representation where we treat p{t) as an operator in Hilbert 
space. In the second equation we consider p{t) as a vector in Liouville space. We further 
introduce the Liouville Space evolution operator 



where LA = [H, A], is the Liouville operator. 

We shall denote superoperators by a subscript u = L,R where the operators Al and An 
act on the ket (left) and bra (right) of the density matrix {A^B = AB and Aj^B = BA)^. 
We further define the equilibrium distribution function 



p(t) = f/(t)MO)f/t(t)^^(t)p(0). 



(21) 




(22) 




(23) 
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Adopting this notation for Eq. (fTTj) yields for the linear response 




(24) 



and for the third order response 



Rrin,Ts,T2,n) = Tr[BLg{m)ARg{T;2)ARg{T2i)ALPeq\ 



(25) 
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Rfin,Ts,T2,n) = TT[BLg{n3)ARgiT32)ALg{T2l)ARp,g] 

4'H^4, ra, T2, ri) = Tl[BLg{n3)ALg{T32)ARg{T2l)ARPeg] 
R?\n,T3,T2,n) = TT[BLg{ns)ALg{r32)ALg{T2i)ALPe,]. 



Note that since the density matrix needs only to be propagated forward, Eqs. \2F^) only contain 
the forward propagator g{t) and not its Hermitian conjugate g^{t), which describes backward 
propagation. This is in contrast with the Hilbert space expression (Eq. ((201)) which contains 
both U{t) and W{t). 

Similar to the wavefunction picture, the response functions may be computed by sums 
over states or by semiclassical wavepackets 

pf{t) = dpdqdp'dc^ I p'q')(p'q' I p^;\t) I pq)(pq | . (26) 

Each term (Liouville space path) in Eq. (j25|l can be recast in the formi^^^ 

Rf{T,,T,,T2.T,) = i:T[Bj^\t)] (27) 

where Pj"""* (t) is a density matrix generating function for path j which can be computed using 
two forward moving trajectories representing the simultaneous evolution of the ket and the 
br a^^i^? . In the wavefunction representation we act on the ket only. Propagating the bra 
forward is equivalent to propagating the ket backward. By keeping track of both bra and 
ket simultaneously we can enjoy the physically appealing all-forward evolution. Since the 
various Liouville space pathways are complex quantities, they interfere when added. This 
interference may result in dramatic effects. 

A systematic approach for computing the response functions will be developed in the 
next section. 

III. SUPEROPERATOR ALGEBRA AND THE TIME ORDERED PERTURBA- 
TIVE EXPANSION OF RESPONSE FUNCTIONS 

In Eqs. (I24|l and ()25|1 we introduced the indices L and R to denote the action of a 
superoperator from the left or the right. In the following manipulations, in particular for 
the sake of developing a semiclassical picture, it will be useful to define their symmetric 
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+) and antisymmetric (z/ 



combinations^^ 



A_=Al-Ar; A+ = - {Al + An) . (28) 

Recasting these definitions in Hilbert space using ordinary operators we get A+X = ^{AX + 
XA) ; A^X = AX — XA, X being an arbitrary operator. Hereafter we sliall use Greek 
indices to denote superoperators A^ witli eitlier u = L, R or u = +, 

Hereafter we shall consider operators that depend parametrically on time. This time 
dependence can be either in the Heisenberg picture A^{t) (Eq. (jHUj) ) or in the interaction 
picture Ai,{t) (Eq. P0|) ). By introducing a time ordering operator T for superoperators in 
Liouville space, we can freely commute various operators without worrying about commu- 
tations. T takes any product of superoperators and reorders them in ascending times from 
right to left. More precisely we define 

A^{Ti)B^{t2) T2 < Ti 

B^iT2)A,{n) n < r2 (29) 

_ l[A,{n)B^{n) + B^{n)A,{n)] = n 

where Ai,{t) is either ^(/(t) or A^{t). T orders all superoperators such that time decreases 
from left to right: The latest operator appears in the far left and so forth. This is the 
natural time ordering which follows chronologically the various interactions with the density 
matrix—. The precise order in which superoperators appear next to a T operator is immate- 
rial since at the end the order will be fixed anyhow by T. For example, T before an exponent 
means that each term in the Taylor expansion of this exponent should be time-ordered. 

We next introduce the Heisenberg picture for superoperators whose time evolution is 
governed by the Liouville operator 

A,{r) = g\r,0)A,g{T,0) (30) 

with 



(31) 



g{T2,Ti) = 6(^2 - n)exp 
Eq. (jnUj) is the Liouville space analogue of Eq. (fTTj). The expectation value of B 

S{t)=TT[Bp{t)], (32) 

may now be represented in a form 



S{t) = (TS+(t)exp 



dr^(r)i_(r) 



(33) 
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The operator B^{t) corresponds to the observation time, whereas A_{Tj) represent various 
interactions with the external field at time r, and (■ ■ ■) denotes averaging with respect to 
the equilibrium density matrix p^q. 

{F) ^ Tr[Fpe,] (34) 

By expanding the exponent in the r.h.s. of Eq. ()33p in powers of E{t), we obtain for the 
response functions. 

i?(")(r„+i ...n)= Q)" (fi+(r„+i)i_(r„) . . . A^n)). (35) 

Eq. ()35j) is merely a compact notation for Eq. ()17|) . It should be emphasized that all time 
arguments are fully ordered ti < T2 . . . < t„_(_i. The Liouville space correlation function in 
the r.h.s. represents a combination of ordinary (Hilbert space) correlation functions. 

Eq. (jH^ may be evaluated directly only for simple models. To convert it into a general 
computational tool we need to develop a perturbation theory for response functions based 
on time ordered superoperators. To that end we partition the Hamiltonian into a simple 
solvable, (usually quadratic) part Hq and a perturbation V 

H = Ho + V, (36) 

and introduce the Heisenberg and interaction pictures. We define the Liouville operators 
L = Lq + corresponding to Eq. (jHUj) where Lq = (Hq)^ i.e., LqX = HqX — XHq. The 
time evolution operator with respect to Lq is 



^o(t"2,Ti) = ^(7-2 - n) exp 



-^^0(^2 - n] 



(37) 



The total (Heisenberg) time evolution operator with respect to L will be denoted G{t2,ti). 
We can then write 

^(T-2,n) = Qoir2,Ti)g{T2,Ti) (38) 

where Q is the time evolution operator in the interaction picture 



G{r2,ri) = Texp 



2 r'^ 



rdT2V_{T) . (39) 
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Throughout this paper we use a hat ( " ) to denote operators in the Heisenberg picture 
(Eq. ipmi ) and a tilde (~) for operators in the interaction picture, i.e. 

A,{r) = gl{T,0)A,goiT,0) (40) 

u = +,—, or u = L, R. 

The equihbrium density matrix of the interacting system can be generated from the 
density matrix of the noninteracting system (po) by an adiabatic switching of the couphng 
V, resulting in 

Peg = GiO, -Oo)po- (41) 

For isolated system at zero temperature, Eq. PT|) generates the ground state density 
matrix of the interacting system, starting with the noninteracting ground state. This is the 
procedure of Gell-mann and Low'^'^. At zero temperature the zero order ground state evolves 
into the actual normalized ground state and hence Eq. ()4H) need not have a denominator. 
Note that in the wavefunction (Gell-mann-Low) formulation of adiabatic switching, the 
wavefunction acquires a singular phase which must be cancelled by a denominator given by 
the closed loop S matrix; the Liouville space expression is simpler since the phase never 
shows up. A remarkable point is that Eq. ()41|) holds as well at finite temperatures provided 
the system is coupled to a bath at constant temperature. This is a thermodynamic adiabatic 
switching where the populations of adiabatic states change and equilibrate with the bath at 
all times^i^^^. It is distinct from the adiabatic switching of an isolated quantum system 
where the populations of adiabatic states do not change^. 

At finite temperatures we start with a grand canonical distribution 

^ exp[-P{Ho-fiN)] 
^° Trexp[-l3{Ho- pN)] ^ ' 

Where (3 = (ksT)^^ {ks is a Boltzmann constant); /i is a chemical potential, N is the 
number operator of particles, and Eq. (PT|) generates the distribution. 

exp[-(3{H - f,N)] 
Trexp[-/3{H-fiN)] ^ ^ 

We now have all the ingredients required for computing the response. Let us start with the 
linear response function 

R^'\T,,n) = '-{B^{r,)A^{n)). (44) 
12 



Using Eqs. and (jHHj) we obtain 

R^'\r2,n) = ^rr[^t(^^^o)fi+(r2)^(r2,0)^t(^^^0)i_(ri)^(ri, 0)^(0, -oo)po] (45) 
The last ^"^(r2,0) can be neglected since it does not affect the trace. Also 



(46) 



which gives 



R^'\r2,n) = -Tr[5+(r2)^(r2,ri)i_(ri)^(ri,-oo)po]. 



(47) 



The time ordering operator allows us to express Eq. ()47|) in the compact form 
R^'\r2,n) = ^^T5+(r2)i_(ri)exp[-^£^rfrV_(r 



(48) 



where we define averaging with respect to the density matrix po of the noninteracting system 



(F)o = Tr[Fpo]. 



(49) 



Eq. ()48|) can be immediately generalized for the response to arbitrary order 



S{t) = {TB+{t) exp 



hJ- 



rfrVl(r) 



exp 



drE(r)i_(r) 



(50) 



Expanding Eq. (fKn|l to n'th order in the external field gives 



R^^\Tn+i---n) = i^-j {TB4Tn+i)A.{Tn)---A.{T,)exp[--J_^ dr\/_(r)])o, (51) 



where we recall that 



X(t) 



exp \ -^Lorj X exp (^-^Lqt 



X = A+,A_,V- 



The Taylor expansion of the exponent in the r.h.s of Eq. ()5H) finally gives 



(52) 



Tn+l, ■ ■ -nj 



m=0 



/Tn + l f'n + l 

-oo J —oo 



X (T5+(r„+0^-(rn) . . .i-(n)V^-(r:) • • • V.{t[))o. 



(53) 



Eq. (|^ constitutes the interaction-picture representation of the correlation function 
Eq. ^4j4i superoperators in this expression should be time ordered chronologically 
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from right (early time) to left (late time). This forms the basis for formulating a field theory 
and Wick theorem in Liouville space in the next section. 

Note that simple time ordering in Liouville space is a more complex operation when recast 
in Hilbert space. This is the essence of why using superoperators makes the bookkeeping 
straightforward. To demonstrate that, let us take -R*-^-* (we use the Heisenberg picture but 
the argument holds as well in the interaction picture, where we should simply replace all 
(")by(~)). 

R^'\T,,T„n) = (0 (TB+(r3)i_(r2)i-(n))o. (54) 

We need to apply the superoperators in a time ordered fashion (in Liouville space) i.e., first 
A_{ti), then A_{t2) and finally B^It^). Separating all possible actions for the left and the 
right we get 

2Tr[T5+(r3)i_(r2)i_(ri)pe,] = (55) 

TT[B{T,)A{T2)A{n)p,g] + TT[A{T2)A{n)p,gB{n)] 

-TT[B{rs)AiT2)peAn)] - Tr[i(r2)pe,A(ri)5(r3)] 

-TT[B{n)A{n)p,,A{T2)] - 'Yx\A[t^)p,,A[t2)b[t-^\ 

+Tr[fi(r3)p,,i(ri)i(r2)] + TY\p,,A{r^)A{T2)B{r^)\ 

In Hilbert space (r.h.s of Eq. (jKK|l ) all operators which act on p^q from the left are time 
ordered and the time increases as we go to the left starting with p^q. All right operators are 
ordered in the opposite way: Time increases as we go to the right starting with peg- This 
mixture of positive and negative time ordering coming from the evolution of the ket (left) 
and the bra (right), respectively is what complicates the bookkeeping of ordinary operators 
in Hilbert space. This is in marked contrast with Liouville space (l.h.s. of Eq. (j55|) ) where 
we keep track of the left and right labels of the various interactions. Consequently all 
superoperators are always positively time ordered in real, physical time which makes the 
formulation of a Wick theorem possible. 



IV. THE CUMULANT EXPANSION AND WICK'S THEOREM FOR BOSON SU- 
PEROPERATORS 

So far we considered four types of operators which enter Eq. ()5H): the reference Hamil- 
tonian — pN\ A, representing the coupling to the external field; V, representing the part 
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of the Hamiltonian to be treated perturbatively, and the desired observable B. To proceed 
further we need to introduce the concept of elementary operators. Any dynamical system 
can ultimately be described by a basic set of operators whose commutators (or anticommu- 
tators) are c numbers. Examples for elementary operators with commuting algebra (EGA) 
are the canonical variables [Qa,Pi3] = i^^ap and Boson operators, [aQ,,a^] = Sap used to 
describe systems of identical Bosons in second quantization. Second quantized fermions are 
described by elementary operators with anticommuting algebra (EAA) {cq, cjj} = 6ai3- The 
operators X = A, B,V, Hq and are some functions of these elementary operators. 

We choose our reference to be a quadratic Hamiltonian given by a bilinear combination 
of elementary field operators. 

Ho = J dxT{x)^\x)ij{x) (56) 
or using creation/annihilation operators 

Ho = Yl Trsolas (57) 

where 

^{^) = ^Vs{x)as (58) 

s 

and ips is a single particle basis set. For Bosons, these operators satisfy the commutation 
relations 

[as, al] = Srs (59) 

and 

[tlj{x),ilj\x')]=6{x-x'). (60) 

For Fermions, Eq. (j59|) should be replaced by an anticommutator. Our elementary set of 
operators is thus the set as,aj or the field operators '?/'(x), ^/'^(x). The following arguments 
hold for Fermions as well, however the derivation is simpler for Boson fields with EGA. We 
shall therefore focus on Bosons first, and the extension to Fermion fields will be presented 
in Section VI. 

We will denote the elementary operators as Qj and introduce the corresponding super- 
operators Qji, V = L,R,+,—. We first note that the superoperator corresponding to any 
function of Qj can be expressed in terms of Qj+ and Qj-, i.e., 

ifiQj)]- = - fiQ.R) = f (g,+ + Iq,-) - f (g,+ - Iq,-) , (61) 
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and 



2 [/(g,)]+ = /(g,x) + fiQ.R) = f [ g,+ + ^Qj- ) + / ( - ^g,- ) . (62) 



(Q?)+ = Q?+ + (63) 



For example 

1 
4 

and 

(g,2)_ = g^.^g._ + g^._g^.^. (64) 

Using these rules (and additional useful relations given in Appendix A) we can expand 
B^{t), A_{t) and V^{t) in a Taylor series in Qj^ and Qj^, converting the time ordered 
product of superoperators in Eq. (|^ into a time ordered product of elementary operators. 
We thus need to calculate 

W{j 

{tn) ■ ■ ■QnuA'^i))o^ (65) 

where I'l,..., i/^r = ± and jm runs over the various operators. The number of operators 
in such products that enter the computation of i?*^"-* is greater than n + 1,N > n + 1. 
The reasons are (i) A^, By may be nonlinear functions of elementary operators and we use 
Eq. (jF)T|l and the formulas of Appendix A to express them as products of Qy. (ii) The 
expansion in V- adds more operators to the product. 

To compute W we define a superoperator generating functional 



5({J(t)})= Texp 



H J Jjuir)Qjy{T)dT 



(66) 



Time ordered correlation functions of superoperators can be obtained from the generating 
functional by functional derivatives 



W{j^U^T^} = ^ . . . -—^—-S{Jit)} 



(67) 

J=0 



Since the Hamiltonian is quadratic, the generating functional may be computed exactly 
using the second order cumulant expansion. This gives 

SiUm) = exp I 5] r dT2 r dn (68) 

[ j,k ''^^ -^-^ 
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where we have introduced the two fundamental Liouville space Green functions. 

G%-iT, - n) = ^{TQ,4T,)Q,_{n))o (69) 

Using Eq. p8|l we can recast these Green functions in Hilbert space 

GtjTir) = '-e{r) [(4(r)Q,(0))o - (Q,(0)Q.(r))o] (70) 

G%^{r) = lm{r)Q,{0))o + (4(0)Q,(r))o]. (71) 

The factor in G^ was introduced for making the classical limit more transparent (see 
next section); since with this factor G~^ has a well defined classical limit. Note that since 
the trace of a commutator vanishes identically, in a time ordered product the superoperator 
to the far left must be a "+". The Green functions G and G ^ thus vanish identically 
and we only have two fundamental Green functions and G^ . Note also that G^^{t) is 
finite for positive and negative r whereas G^^(r) vanishes for r < 0. Eq ()68|1 is an extremely 
compact expression which can readily be used to compute response functions to arbitrary 
order. 

The two fundamental Green functions can be expressed in terms of the matrix of spectral 
densities Gij{uj) defined as the Fourier transform of (7+ -^3^4i^52^67 . 

—G.,{u)sin{ur). (72) 

-oo ZTT 

We then have 

The Wick theorem for superoperators then follows from Eqs. (j67|) and (|68|) and can now be 
stated as follows: 

{TQnuAn) . . . Qi^^^(rjv))o = (74) 

V 

Here jai^a . . .jqVq is a permutation of jiz/i . . .j^v^ and the sum runs over all possible per- 
mutations, keeping the time ordering. Since only G^^ and G^~ survive many of the terms 
will vanish. 
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Wick's theorem makes it possible to develop Feynman diagram perturbative techniques 
which express the linear and non-linear responses of the interacting system in terms of the 
two fundamental Green functions. This theorem is useful whenever a quadratic reference is 
adequate and non quadratic parts of the Hamiltonian can be treated perturbatively. It states 
that multipoint correlation functions of systems with quadratic Boson Hamiltonians may be 
factorized into products of two-point correlation functions of the primary coordinates. 

V. MODE COUPLING AND SEMICLASSICAL RESPONSE OF BOSON FIELDS 

Eq. ()17|) contains 2" terms representing all possible "left" and "right" actions of the 
various commutators. Each term corresponds to a Liouville-space path and can be repre- 
sented by a double-sided Feynman diagram^S,. The various correlation functions interfere 
and this gives rise to many interesting effects such as new resonances. The {i/h)"' factor 
indicates that individual correlation functions do not have an obvious classical limit. The 
entire response function must, however, have a classical limit. When the various correlation 
functions are combined, the (i/h)"' factor is cancelled as h tends to 0, and one obtains the 
classical response, independent of h. The elimination of h for higher nonlinearities requires 
a delicate interference among all 2" correlation functions. 

The terms contributing to -R*-"-* (Eq. will generally have a {i/fi)'^^^ factor where p 

is the order in V-. This factor must be cancelled as ^ to ensure a well defined classical 
limit. This is guaranteed since by Wick's theorem we will have n + p terms, each 
carrying an h factor. In the classical limit we set cothlhu /2k bT) = 2kBT/huj. We then see 
from (Eq. [721) ^^at the two Green functions are simply connected by the classical fluctuation 
relation. 



G^^ is independent on h. The factor hcoth.{huj/2kBT) = h/ taiih.{huj/2kBT) is analytic 
in h and can be expanded in a Taylor series, thus yielding a semiclassical expansion of the 
response. To obtain the classical limit we need to keep h in the generating functional, perform 
the h expansion (since response functions are generally analytic in h) and only then send 
h —>■ 0. Setting this limit at the right step is essential for developing a proper semiclassical 
expansion. Classical response functions may not be computed using classical trajectories 
alone: The response depends on the vicinity of a trajectory. One needs to run a few closely 




(75) 
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lying trajectories that interfere. Formally this can be recast using stability matrices which 
carry the information on how a perturbation of a trajectory at a given time effects it at 
a later time. The repeated computation of the stability matrix greatly complicates purely 
classical simulations^i^. The semiclassical expansion circumvents these calculations in a 
very profound way. Corrections to the trajectory to low order in h carry the necessary 
information. Combining several semiclassical trajectories^ allows them to interfere and the 
leading order in h for i?") survives and gives the classical response. This allows us to 
avoid computing stability matrices which is required when the classical limit is considered 
from the outset. The classical limit obtained this way reproduces the results of mode coupling 
theory and removes all ambiguities as to how higher order correlation functions should be 
factorizecii^'i^'i^. 

To illustrate how this works let us consider the following model Hamiltonian Hm 

25.26.41.51.52 



/^..-E^^ + Mj+V-IQ). (76) 



where Pj{Qj) is the momentum (coordinate) operator of the j'th primary mode, Qj and Mj 
are its frequency and reduced mass respectively and ^(Q) is the anharmonic part of the 
potential. The primary modes interact with a large number of harmonic (bath) coordinates 
which induce relaxation and dephasing. Low-frequency bath degrees of freedom q and their 
coupling to the primary modes are described by the Hamiltonian Hb and the material 
Hamiltonian is given by22*2^ 

H = HM)+Hb{Q,ci). (77) 



We assume a harmonic bath linearly coupled to the primary coordinates Qj 



2m„ 2 \ "^o^a 



(78) 



where PaiQa) are momentum (coordinate) operators of bath oscillators. 

This model gives the following Brownian oscillator form for the spectral density^ 

M, VL and E are diagonal matrices and their matrix elements are Mij = SijMj, Qij = SijQj 
and lij = 6ij. C"{u!) is the odd part of the spectral density, which is related to the even part 
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C'ioj) by the fluctuation-dissipation theorem. The spectral density (Eqs. (f7^ and (f75j) ) is 
given by C{u) = [1 + coth(/?^a;/2)]C"(cj). 

(Sjj) is the imaginary (real) parts of a self energy operator representing relaxation 
(level shift). 

T.,Auj) = --Re r du'^^^, (80) 
TT J-oo cj' — u; 

Equations and (j^ . together with ((721), and (f?^ constitute closed expressions 
for the Brownian oscillator response functions. Ordinary Langevin equations are obtained 
by taking the overdamped limit 7 >> Vt oi Eq (f7^ . When the primary coordinates are 
uncorrelated the superoperator Greens functions are 

Gtf{r) = 20(r) exp(-A,r)A,Ai<5,, (81) 
G++(r) = ;iexp(-A,r)A,A, coth(/3^A,)5,,- + 1 £ ''-IL^^^t^X.-XS,, (82) 

P n=l '^n ^^i 

where Vn = 27m /hf3 are the Matsubara frequencies, A, = Qf/'yu and Aj = l/Mjfi?. The 
expansion of nonlinear response functions using collective coordinates have been discussed 
in detai P^i^^i^^i^'^ and recently employed in mode coupling theoryi^*i^. 

All nonlinear response functions of the linearly driven harmonic oscillator vanish iden- 
tically due to interference among Liouville space paths-^. The simplest model that shows 
a flnite nonlinear response is a nonlinearly driven Harmonic oscillator where the operator 
y4 is a nonlinear function of the coordinate. This model has been studied both quantum 
mechanically and classically^. Its response can be alternatively computed by following the 
dynamics of the Gaussian wavepackets in the complete (system and bath) phase space, since 
the system-bath Hamiltonian Hb is harmonic in the full phase space {P, Q, pa, Qa] ^^i^^i^^i^^i^^ . 

We next discuss the connection between our results and a fully classical computation of 
the response. In classical mechanics the density matrix assumes the form of an ordinary 
distribution function in phase space. This can be obtained from the quantum density matrix 
by switching to the Wigner representatioii^. 

Pw{p<i;t) = j^^j^ J rfs(q- s/2|p(t)|q + s/2)exp(p.s) (83) 

where d is the number of degrees of freedom. The Wigner representation offers a transparent 
and simple semiclassical picture that interpolates between the quantum and classical regimes. 
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Wavefunctions, on the other hand, do not have a clear classical counterpart (although there 
are, of course, very powerful semiclassical approximations for the wavefunction such as the 
WKB approximation). 

Wick's theorem for superoperators in Liouville space allows to develop a unified picture 
of quantum field and classical mode coupling theories which clearly reveals the information 
content of classical and quantum nonlinear response functions. Both classical and quantum 
response functions contain interference. Quantum mechanically it is between 2" Liouville 
space paths. The classical interference is of a very different nature^and involves 2" closely 
lying trajectories. The response function in phase space is obtained by an ensemble averaging 
over such bundles of trajectories^!^-. Alternatively, the classical response can be recast using 
stability matrices which carry the relevant dynamical information on the vicinity of a given 
trajectory. The connection between the quantum and classical 2'^-fold interference is more 
transparent by keeping the left /right or the +/— pathways rather than working in phase 
space. We keep % alive during the semiclassical calculation and send it to zero only at the 
end. 

In the fully classical phase space approach, we take the two separate (left and right) paths 
required in a quantum mechanical formalism and expand them around a single classical 
reference path, letting the stability matrices carry the burden of retaining the information 
about the differences between the paths. In the present +/— representation we keep track 
of closely lying trajectories by retaining terms to order h and combine them only at the very 
end. This way stability matrices which pose enormous computational difficulties^ never 
show up. 

Another way to view the classical/quantum connection is by starting with the expressions 
for quantum mechanical nonlinear response functions in terms of combinations of n point 
correlation functions of the relevant variables. These correlation functions differ by their 
time ordering i.e. {A{ti) A{t2) A{T-i)) ^ {A{t2)A{ti)A{t^)) etc. i?*^"^ is then a combination of 2" 
such correlation functions, each representing a distinct Liouville space pathway. Classically, 
of course, time ordering is immaterial since all operators commute and we have only a 
single n point correlation function. The presence (absence) of symmetry with respect to the 
permutation of the n time variables in classical (quantum) correlation functions implies that 
the effective multidimensional space of time is reduced by a factor of n\ in the classical case. 
Classical correlation functions thus carry less information than their quantum counterparts. 
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Classically, it suffices to calculate {A{ti)A{t2) ■ ■ ■ A{Tn)) for ri < r2 ■ ■ ■ < t„. Quantum 
mechanically, each of the n\ permutations of the time arguments is distinct and carries 
additional information. The stability matrices provide the extra information required for 
computing the response functions from classical correlation functions. 

Since classical correlation functions do not carry enough information for computing non- 
linear response functions, it is not possible to simulate and interpret the response in terms 
of standard equilibrium fluctuations; additional nonequilibrium information^ is necessary. 

Correlation functions are equilibrium objects and can be computed using sums over un- 
perturbed trajectories; response functions can be obtained either as equilibrium averages 
(stability matrices) or recast in terms of 2" closely lying nonequilibrium trajectories per- 
turbed at various points in time. Quantum corrections to classical response functions may 
be represented in terms of classical response functions of higher orders^. 

Finally we note that an alternative semiclassical h expansion of the response is possible 
even when the temperature is low compared with the material frequencies, and the system is 
highly quantum, provided anharmonicities are weak— >2i. The leading terms in the expansion 
can be obtained by solving classical equations of motion. This is done by hiding the h in the 
coth factor in Eq. (f72|) by recasting it in the form coth(ti;/2ci;r) where ujt = ksT /h is the 
thermal frequency, and redefining G^^ by multiplying it by h. The response is then analytic 
in h (as long as we forget about the h dependence of ojt)- Semiclassical approximations 
ordinarily hold when the temperature is high compared to all relevant vibrational frequencies. 
This points to a much less obvious, low-temperature weak anharmonicity regime, where the 
response of the system behaves almost classically even though its temperature is very low. 



VI. WICK'S THEOREM FOR FERMION SUPEROPERATORS 

One reason why the handling of Boson operators and fields is simpler compared to 
Fermions is that superoperators corresponding to elementary Boson operators are also ele- 
mentary i.e., their commutators are also numbers. To see that let us consider the commu- 
tation rules of superoperators by acting with commutators on an arbitrary operator F. 

[Q,L,QkL]F = [Qj,Qk]F (84) 
[QjR-,QkR\F = —F[Qj,Qk] 
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[QjL, QkR\ — 

Since [Qj, Qk] is a number, we see that the corresponding superoperators are elementary as 
well. This property holds also if we consider the linear combinations in the +/— represen- 
tation. To see that we start with [y4+, BJ\ and act on an ordinary operator F 

= \[A,B]F+\f[A,B] (85) 
[A.,BJ\F = 4[A+,i?+]F= 

Since the commutator of elementary operators is a number, Eq. (j85|) gives 

= [A,B] (86) 
= 4[A+,fi+] = 

The commutators of elementary Boson superoperators are thus numbers. It then follows 
that the superoperators corresponding to elementary Boson operators are Gaussian, be it 
in the +/— or the L, R representation. The indices v in Eq. (|7^ can thus run over +/— or 
L, R and Wick's theorem holds in either case. 

Using L, i?, the functional can be used to generate individual Liouville pathways. Using 
H — it generates combinations of such paths, making the classical limit more transparent, 
since we work with combinations of Hilbert space correlation functions which enter the 
response and have well defined classical limits. 

Life is more complicated for Fermions. To see that let us consider the anti-commutation 
rules for Fermi elementary superoperators. 

{QjL,QkL]F = {Qj,Qk]F (87) 
{QjR,QkR}F = F{Qj,Qk} 

{QjL, QkR}F = 2QjFQk 

Since {Qj,Qk} is a number this shows that left/right or right/left superoperators have ele- 
mentary anticommutators but the anticommutator of left / right operators is not generally a 
number. We thus do not have a Gaussian statistics. Note however that left and right oper- 
ators always commute. The generating functional needs to be formulated using Grassman 
algebra of anticommuting numbers similar to standard Green functionsi*^. The important 
point is that a modified Wick theorem still holds for Fermi superoperators. It is given by 
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Eq. ()74|) with the following changes: 

(i) We must use L, R rather than +/— variables for u. 

(ii) Each term needs to be multiplied by (—1)^ where P is the number of permutations 
of elementary operators required to bring them to the specified order. Since left and right 
superoperators commute, we only count the number of permutations among left and among 
right operators. (Permuting a right and left operator does not count in P). 

The expectation of the T product of any number of (Boson or Fermion) superoperators 
may thus be expressed as the sum of all possible products of expectations of T products of the 
separated pairs of operators for the reference many-particle density matrix po corresponding 
to Ho. 

VII. DISCUSSION 

Hilbert space and Liouville space offer a very different language for the description of 
nonlinear response. The density matrix provides a fully time ordered description since we 
only need to propagate it forward in time. In contrast, the wavefunction involves both 
forward and backward propagations. The choice is between following the ket only, moving 
forward and backward or following the joint dynamics of the ket and the bra and moving only 
forward. Artificial time variables (Keldysh loops) commonly used in many-body theory22iSL^ 
are connected with the wavefunction picture. The density matrix uses the real laboratory 
timescale throughout the calculation. 

In Liouville space all observables are time-ordered leading naturally to a semiclassical ap- 
proximations^- and Feynman path integral diagrammatic techniquesii^iMi^iSil. Maintaining 
time ordering allows to recast S**-"-* using nonlinear response functions which decouple the 
field and the material parts. The nonlinear response function is calculated as a path integral 
in Liouville space by summing over the various possible pathways in Liouville-space which 
contribute to the polarization. Path integrals have been extensively used as a useful tool 
for numerical computations of mixed quantum-classical calculations^i^. The density ma- 
trix formulation provides a similar development for phase space-based numerical procedures. 
Graphic visualization of these paths is provided by double-sided Feynman diagrams^^. 

The density matrix Liouville Space picture offers many attractive features. Physical 
observables are directly and linearly related to the density matrix. Consequently, every step 
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and intermediate quantity appearing in the description has a simple physical meaning and 
a clear classical analogue. This should be contrasted with wavefunction-based calculations 
of a transition amplitude, which by itself is not an observable since signals are related to 
sums of products of such amplitudes. 

The density matrix provides a practical way for performing ensemble averagings and de- 
veloping reduced descriptions where bath degrees of freedom are traced out from the outset. 
Since it represents the state of the system by a matrix rather than a vector, an N point 
grid for p and q in a semiclassical picture will require points for the wavefunction and 
N"^ for the density matrix. The ability to perform ensemble averagings and reduced de- 
scriptions more than compensates for this price for complex systems. Many-body theory 
of superoperators further naturally allows for the treatment of dcphasing and decoherence 
effects. Diagonal and off diagonal elements of the density matrix are known as populations 
and coherences, respectively. When an off diagonal element evolves in time for a system 
coupled to a bath, it acquires a phase since its evolution from the left (ket) and the right 
(bra) is governed by different bath Hamiltonians. This phase depends on the state of the 
bath. When we perform an ensemble average of these elements over the distribution of 
the bath degrees of freedom, this variable phase results in a damping of these elements. 
The damping of off diagonal elements of the density matrix resulting from phase (as op- 
posed to amplitude) fluctuations is called pure dephasing or phase relaxation (also known 
as decoherence). Dephasing processes can only be visualized in Liouville space by following 
simultaneously the evolution of the bra and of the ket and maintaining the bookkeeping of 
their joint state. Dephasing processes directly affect all spectroscopic observables since they 
control the coherence which is the window through which the system is observed. Differ- 
ent pathways representing distinct sequences of populations and coherences, are naturally 
separated in Liouville space. 
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APPENDIX A: SOME USEFUL RELATIONS FOR SUPEROPERATOR ALGE- 
BRA IN LIOUVILLE SPACE 



A Liouville space operator is labelled by a Greek subscript where a = L,R, +, — . It is 
defined by its action on X, an ordinary (Hilbert space) operator. We write a general matrix 
element 

{A^X)ij = J^iAahM^Ki (Al) 

Aa is thus a tetradic operator with four indices. Since AiX = AX and AjiX = XA we 
obtain using Eq. (|A1|) 

{AL)ij,Ki = Ai^5ji (A2) 

{AR)ij^^i = Aij5if, (A3) 

Note that the order of the ji indices in Eq. flA3|l has been reversed. 

Since A+ = ^{Al + Ar) and A_ = Al - Ar we have (A_)jj-^£ = AiJji - Aij6i^ and 
{A+)ij^f^i = ^[Aif,6je + A^jSi^]. It then follows that [Al,Br] = 0. This commutativity of 
left and right operators is possible thanks to the large size of Liouville space and simplifies 
algebraic manipulations resulting in many useful relations. 

2[A+, 5_] = [Al, Bl] - [Ar, Br] = [AB)^ - [BA)^, (A4) 
exp(Ai) = exp(yl+)exp(|y4_), 
exp(A+) = exp(iAL)exp(|Afl), 
exp(yl_) = exp(yli) exp(-AR), 
(expA)+ = 2exp(A+)cosh(iA_), 
(expy4)_ = 2exp(A+) sinh(|y4__). 

In the following a is a complex number 

6{cu - A_) = j da 5{a- Al)5{ijJ - a + Ar), (A5) 

5{ijj — A_|_) = j da 5{a — Al)5{(jJ — a — Ar), 
5{a - AL)5{a - Ar) = 6{A+ - a)6{A_). 
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